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Abstract 

In this note we present a series expansion of inverse moments of a non-negative discrete random variate 
in terms of its factorial cumulants, based on the Poisson-Charlier expansion of a discrete distribution. 
We apply the general method to the positive binomial distribution and obtain a convergent series for its 
inverse moments with an error residual that is uniformly bounded on the entire interval < p < 1. 
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1 Introduction 

Certain problems in statistics and in other branches of science require the calculation of the inverse moments 
of a distribution, also called the negative or reciprocal moments. In statistics, noteworthy examples are life 
testing problems [10] and sampling problems with samples of random length |22| . Especially useful are the 
inverse moments of the binomial distribution. Recent applications in quantum physics, for example, include 
the calculation of running times of certain quantum computation algorithms |33| , the study of random walks 
on n-dimensional cubes and exact calculations of the confidence region of a Beta estimator of the 
parameter p of a binomial distribution [2] • 

While inverse moments of positive binomial variates can be calculated exactly [28j , the ensuing expressions 
are (A^-|- l)-term summations and the calculations become complicated for large TV. In addition, it is far from 
clear from these summation formulas what the asymptotic behaviour with N should be. For these reasons 
there is an interest in obtaining efficient series expansions. Inverse moments of positive binomial variates 
have been studied as early as 1945 by Stephan [28], who considered expected value and variance of negative 
powers and expressed them as series expansions of inverse factorials. In that work an even earlier reference 
was made to the work of Bohlmann [6l , whose approach was to expand the function 1/a; in a Taylor series 
and take expected values of each term. This approach has recently been revived in Ref. [32] . Based on these 
expansions, Grab and Savage [13] have calculated tables for the first inverse moment of positive binomial 
and Poisson variates. Many other expansions have been proposed, for example expansions in Eulerian 
polynomials |22j and in factorial powers [25j . These expansions also work for related distributions like 
the inverse binomial and Poisson distributions. Govindajarulu [12] has found recurrence relations between 
inverse moments of positive binomial variates and Refs. [T71 [HI [301121112? contain various bounds on inverse 
moments. More general methods, valid for any distribution and involving integrals, have been proposed 
in [71 [11 [9l [181 [26] . The problem with any of these series expansions considered so far is that as concerns 
convergence they do not perform equally well over the complete range < p < I. Some of these expansions 
are asymptotic series and are divergent. We illustrate this claim with graphical examples in Section [3] 

In this paper we present a new series expansion for inverse moments, based on a very simple idea. A 
discrete distribution of a non-negative discrete random variable can be approximated by a series expansion 
based on the Poisson distribution [3] . The coefficients of this expansion arc the factorial cumulants of the 
original distribution. This expansion then induces an expansion of the inverse moments in terms of the 
inverse moments of the Poisson distribution. The main point we wish to make in this paper is that this 
expansion of inverse moments is a very good one because of its excellent convergence properties, converging 
very rapidly and consistently throughout the interval < p < 1. 

In Section [4] we give a brief overview of the Poisson expansion method, and in Section [5] we derive 
the required expectation values of the Poisson distribution in terms of its inverse moments, which can be 
calculated with standard software or alternatively via the formulas presented in Appendix [X] Our main 



result is Theorem 1, Section [51 The excehent convergence properties are obvious from Figures [3] and [51 
which show the absolute and relative error of the expansion applied to the binomial distribution. 



2 Notations 



Throughout this paper, Q denotes a Poissonian random variate with mean value fi. Its probability dis- 
tribution function (PDF) is given by the semi-infinite sequence tt^ = (7r^(0), 7r^(l), . . . , 7r^(fc), . . .) with 
Tr^{k) = e-^///fc!, for fc = 0, 1, . . . 



Let / be a semi-infinite sequence 



/-(/(0),/(l),...,/(fc),...). 



We will always set /(fc) = for A: < 0. The forward difference operator A and the backward difference 
operator V are defined via 

Af{k) = /(fc + i)-/(fc) 
Vf{k) = f{k)~f{k~l). 

These operators can be represented by semi-infinite matrices: 



-1 1 

-1 1 



V V 



f 1 ^ 

-1 1 

-1 1 



V 



Higher-order difference operator A' and V' {I > 1) are defined as 



I 



J=0 



These operators are represented by the Z-th matrix powers of A and V: 

A' ^ A' 

V' ^ V' = (-l)'(A')^. 



(1) 
(2) 



The factorial cumulants k'^' of a discrete distribution with PDF / are generated by the logarithm of the 
expectation value of (1 -|- x^: 



\fe=0 / j=0 



(3) 



To calculate certain inverse moments of the Poisson distribution, we will need the Stirling numbers of 
the first kind S'j*^'' [T]. They satisfy the recurrence relation 



sf^^ = sf-^^ -jSf\ j>k>l, 



(4) 



with s[^^ — 1, and are generated by 

n — 1 n 
*:=0 j=l 

We also need the non-central Stirling numbers of the first kind Sj'^j^ [TSl [20] ; for Z = they coincide with the 
ordinary Stirling numbers. They satisfy the recurrence 

= j>k>i, (6) 

with 5'!^^^ = 1, and are generated by 

n—l-\-l n 

k=i+i j=i 



For k — I explicit formulas exist: 
and 



sl'^ = {-iy-\j-i)i (8) 
4V = (-ir'(j + ^-i)!A! (9) 



3 Numerical comparison of existing expansions 

Here we consider three existing expansions of the first inverse moment of the binomial distribution and 
illustrate their convergence behaviour over the entire interval < p < 1. As is customary, we put q = I — p. 
The oldest known expansion is Stcphan's expansion "28]. The M-term expansion reads 

il-q-)E[l/K\K>0] = t ^^f^ (10) 

:= l-±(''^')p^,---^. (11) 
j=o V / 

For very small p, this expansion suffers from numerical instability. For moderately small p, convergence is 
very slow, as can be seen from Figure [1] 

A recently obtained expansion is Rempala's |25j , which is an improvement on an expansion by Marciniak 
and Wesolowski [25] : 

(1 - q'^ni/KlK > 0] = (Np)-' ^ (12) 

i=0 \ i ) 

For not too small p this series seems to converge much faster than Stephan's expansion. However, as it is 
an asymptotic expansion, it actually diverges. For fixed p there is an optimal number of terms M , and for 
larger M the error increases very rapidly. For example, taking N ~ 100, the first 100 terms of the series give 
extremely accurate results for p > 0.54 but are completely useless for p < 0.51. The absolute error of this 
expansion is illustrated in Figure [U 

Znidaric [32] gives an expansion formula for all inverse moments. When specialised to the case of the 
first inverse moment it reads: 

(1 - ,"mWK > 01 ^ I' (-!)'(= + Dff^, (13) 



where fii{N — 1) is the i-th central moment of Bin(A^ — l,p)- Just like Rempala's expansion, this is an 
asymptotic series and suffers from the same divergence problems; in addition it gives much less accurate 
results, as seen from Figure [31 
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Figure 1: Absolute error as a function of p of Stephan's expansion of the first inverse moment of a positive 
binomial variate for N = 10 (upper graph) and N = 100 (lower graph), with 1 term (upper curve), and up 
to 10 terms (lowest curve). 
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Figure 2: Absolute error as a function of p of Rempala's expansion of the first inverse moment of a positive 
binomial variate for = 10 (upper graph) and N = 100 (lower graph), with 1 term (upper curve), and up 
to 10 terms (lowest curve). 
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Figure 3: Absolute error as a function of p of Znidaric's expansion of the first inverse moment of a positive 
binomial variate for iV = 10 (upper graph) and N = 100 (lower graph), with 1 term (upper curve), 3,5,7,9 
and 11 terms (lowest curve). Expansions with an even number of terms are left out for clarity. 



4 Poisson expansion of a probability density 



It is well-known that the binomial distribution Bin(Af,p), with probability distribution /(fc) — (^)p'^(l — 
p)""*^, k = 0, . . . , N, tends to the Poisson distribution with mean value /i = Np, when Np is kept fixed and 
N tends to infinity. In [3l [4] Barbour and coworkers presented a series expansion for a general distribution 
of a non-negative discrete random variate in terms of the Poisson distribution, providing a quantitative and 
more generally useful version of this last statement. They gave this expansion the name of Poisson-Charlier 
expansion. 

The first term of the Poisson expansion of a distribution / is, of course, the Poisson distribution itself, 
with mean fi given by the mean of /: 

oo 

/i = E[fc]=^/(fc)fc. (14) 

The higher-order terms of the Poisson expansion consist of backward differences (V) of tt^ (fe) , with coefficients 
based on the factorial cumulants k^''^ of /. Formally, the Poisson expansion is based on the following identity: 



\k=2 



/ = e^p(E V^'^^'l^'^- ^^^^ 



A simple proof of this identity is given in Appendix [B] 

In practice, the right-hand side of has to be replaced by a finite series, with a finite number of terms 
in powers of V, and with suitable error bounds estimating the truncation error. There are many ways to do 
this. The most obvious way is to replace the right-hand side by a Taylor series in V. An m-th order series 
can be obtained by expanding the function 

(°° (k) \ 

in powers of t up to t™~^ around t = 0, and then putting t = 1. This yields a polynomial in V of degree 
TO — 1. We will generally define the order of the expansion as 1 plus the highest degree of t, and its degree 
as the highest degree of V. While conceptually simple, we are not aware of known error estimates for this 
kind of expansion. 

In contrast, Barbour's expansion amounts to an TO-th order series of the function 

i^expj^-^:— HV)'=j, (17) 

which has an extra factor \/t in the argument of the exponential function. Taylor expansion in t yields 
a polynomial of degree 2 (to — 1) in V, hence the degree of the TO-th order expansion is 2 (to — 1). When 
comparing series (|16p and (fT7|) of the same degree one finds that the latter contains a subset of the terms of 
the former. For example, the third order series (I17p is given by 

while the fifth order series fTB]) . which is also of degree 4, has an additional term: 

K-(2) K-(3) /('k-(2)^2 (4)\ 

. . V^- Vv'+(^^ + ^) V^. (19) 



2 6 V 8 24 _ 

Surprisingly, however, while this observation makes the series (|16p look more powerful than (jl7p (more 
terms for the same degree) , the latter is much better suited to the purposes of this paper because it actually 
converges much faster and, moreover, explicit error bounds for it are known. 

The complete expression of (jl7p . including all high-order terms, is quite complicated and will not be 
given here (see [5], eq. (2.7)). Applying this operator polynomial to 7r^(fc) yields the TO-term approximation 
of /(fc). General bounds on the absolute error are given in These bounds were obtained in a highly 
non-trivial way (not to say magical way), using the so-called Stein's method [5]. We will not give the general 
bounds here but will only mention them (in Section [S]) for the special case of the binomial distribution. 



Example 1. The Poisson expansion of the Poisson distribution triviaUy reduces to the first term only, as 
it should, because all factorial cumulants of the Poisson distribution are 0, except k^^^ which is equal to the 
mean value fi. 

Example 2. A more enlightening example is the Poisson expansion of the binomial distribution K ~ 
Bin(iV,p). Its mean is = Np and E[(l + x)^] = (1 + px)'^ , so that 

logE[(l +a;)^] iVlog(l = V (pxY . 

Hence, the j-th factorial cumulant (for j > 0) is 

= -iV(j-l)!(-p)^ (20) 
Thus, the third order Poisson expansion of the PDF of K is 

/(fc) « TT.ik) - ^ vV^(fc) - ^ vV^(fc) + ^ V\,{k) (21) 

with /i — Np. 



5 Inverse moments of the Poisson distribution 

Using the Poisson expansion, calculating the r-th inverse moment of a distribution is reduced to calculating 
the expectations W^l / {Q + aY\Q > Q] for a = 0, 1, . . . where Q is a positive Poisson variate with mean value 
/i. The expectations for a = are essentially the inverse moments of the positive Poisson distribution. The 
other expectations are the inverse moments about —a and can be derived from the moments about the origin 
using a simple recurrence relation. 

To simplify notations we introduce the symbol E+: for any function and for K a positive random 
variate with probability distribution /: 

oo 
k=l 

= (l-/(0))E[5(if)|if>0] 
= E[g(i^)] - /(O)g(O). 

Thus, in particular, 

E+[1/(Q + aY] = E[1/(Q + aY] - e-^^jaF. (22) 
Let's consider first the inverse moments about the origin, E+[l/Q'^]. For r = 1, we have 

oo ^ 

E+[l/Q] = e-^Y.-^ (23) 



gMt _ 1 

At (24) 

The integral is essentially equal to the exponential integral Ei [1]. With 7 w 0.5772 . . . the Euler-Mascheroni 
constant, 

'^^^~T^ = Ei(^)-log/x-7=:Er(/i). (25) 

The exponential integral function is well-studied, and implementations are incorporated in many numerical 
and algebraic software packages. 

The higher inverse moments about the origin can be expressed in terms of hypergeometric functions 



(26) 



which is actually just a standardised restatement of the definition. Again, these functions can be accurately 
and efficiently calculated using standard software packages. In absentia, a second option is to resort to 
explicit series expansions; some of these are given in Appendix [Xl 

Now we move on to the (non-central) inverse moments about —a, E[1/(Q + aY], for a > 0. Let's first 
consider the case r = 1. For a = 1, 

E[f/(g-|-l)] = l^^, (27) 
A* 

and for a > 1 we can use the recurrence [7] 

E[l/(g + a)] = i (1 - (a - l)E[l/(g + a - 1)]) . (28) 
This recurrence can be generalised to higher r as follows: 

Lemma 1. For Q a positive Poisson variate with mean fi, and for integers a > and r > 1: if a = 1, 

E[l/{Q + lY] = -E+[l/Q^-% (29) 

^^ 

and if a > 1, 

E[1/(Q + aY] = - (E[1/(Q + a- lY'^] - (a - 1)E[1/(Q + a - 1)']) ■ (30) 

In [16j a similar recurrence formula is used. 
Proof. First note 

°° .,k 



Ell/(0 + an . «-E^(^ 



fc=0 

" k + 1 



e 



^ (k + 1 

k=0 ^ 



{k + iy. {k + aY 



k + 1 
(fc + 1)! {k + aY 



M t,ik + l) 

^ oo 



oo I: 



fi k\ {k + a-lY 

= -E+[Q/(Q + a-l)'']. 
M 

For a = 1, this directly gives (I29p . 

For a > 1, the fc = term is zero and can be added to the sum: 

k 



= -E[Q/(Q + a- 1)'-]. 



Writing Q = (Q + a - 1) — (a — 1) then yields 



E[l/{Q + aY] = -E[0/(0 + a - 1)'-] 
A* 

= - (E[1/(Q + a- lY'^] -{a- 1)E[1/(Q + a- 1)^) 
/J 



□ 



By iterating the recurrence relations (l27l) . ((28| . (|29|) and pO|) . one can find an expression of the r-th 
inverse moments of the positive Poisson distribution in terms of central and non-central Stirling numbers. 



Proposition 1. Let Q be a positive Poisson variate with mean value /i. For integers a > 1 and r > 1: 

E[1/(Q + aY] = (si^^ (1 - e-n + ^ Sf^^ + ^ si%^, A . (31) 

' V k=l k=l / 

Proof. Let's denote the right hand side of (|3ip by A{a,r), and the r-th inverse moment of the positive 
Poisson distribution E+[l/Q''] by /Lt-r- We consider first the case r ^ a — 1. Since A{1, 1) = (1 — e^'')//i it 
coincides with E[1/(Q + 1)], by ([?f|) . Next, for r = 1 and a > 1, we need to check the recurrence (pS)) . i.e. 
whether 

Aia,l) = {^-{a-l)A{a-l,l))/^l. 
From the generating function ([7]), we have S^^}). j, = (— l)"^'''^^(a — l)!/fc!. Thus, 



^(a,l) 



(-!)-!(« 




It is straightforward to check the recurrence from the latter expression. 

For r > 1 and a = 1 we find A{l,r) = {T,lZ\ ^^''"''V-fc + S^{\l - e"^))//!. As s'f^ ^ 1 if and only if 
r = 1, this simplifies to as required for (|29p. 

Finally, to check the remaining recurrence (|30p for a,r > 1, we first note that the k — a — 1 sum term in 

the last term of ([31]) vanishes, because s[^l = for r > 1. Hence the upper summation limit can be replaced 
by a — 2. Then, from recurrences ^ and (O, we see that the coefficients appearing in the three summation 
terms obey one and the same recurrence: si^~'^^ = sl[~i~^^ ^ {o- ^ ^)Sa-i'\ si^^ = S^~i^ ^ ^ 

and <S'^'2fe k ~ '^a'-i-k fe ^ ('^ ^ l)*^!-!-*; k- then an easy matter to verify that A{a,r) indeed satisfies 
the final recurrence A{a, r) — A{a — 1, r — 1) — (a — l)A{a — 1, r). 

Since A{a,r) satisfies the same recurrences and boundary conditions of E[l/((3 + aY\, the two must 
coincide, proving equality in pT|) . □ 



Inserting the explicit formulas ([8]) and ([9|) immediately gives, for the case r = 1 

(a-l)!(-l^°-i 



[l/iQ + a)] = ' I 1 - + > J-m)Vj! I • (32) 



a-l 

E( 



6 Poisson Expansion of Inverse Moments 

Based on the results of the two previous Sections, we are now in the position to formulate our main result. 
While our own interest lies with the binomial distribution, our result is generally valid for any positive 
discrete random variate. 

Theorem 1. Let K be a positive discrete random variate with probability distribution f , having mean value 
/i and factorial cumulants k*-*-' . 

Let fi-r be the r-th inverse moment = E^"[l/(5''] of a positive Poisson variate Q with the same mean 
value fi as K , and let ^-r,a be the corresponding shifted inverse moments fi-r.a = E[l/((5 + aY], for a > 0, 
obtainable e.g. from \31\} . Define the sequence q-r with q-rio) — fJ'-r,a, for a > 0. 

Let Pm be the degree 2(m — 1) polynomial obtained from the m-th order Taylor approximation of |j7p as 
indicated in Section^ 

Then the m-th order Poisson approximation of the r-th inverse moment of K is given by 

E+[l/i^n « (P„(-A)q_,)(0). (33) 



Proof. Since the Poisson expansion of a distribution is expressed in terms of V^tt^, we first calculate the 
following sums: 



oo 



fc=i 



.a, 

k=l a=0 ^ ^ 

3 



k=l a=l ^ ' k^" 



7r^(fc - a) 



Now note 



E"" 7rp(fc - a) _ 7r^(fc) 

fc'' ~ ^ (fc + aV 

k=a k=0 ^ ' 

= E[l/{Q + aY 



Thus 



3 

J 



E+[l/Qn + ^r (-irE[l/(Q + an 

a=l ^ ^ 



We can express the last equation in terms of the sequence g-,. and the j-th forward difference operator: 

j^^r,, - ((-A)^g_,,)(0). 
Combining this with the m-term Poisson expansion of the distribution of K , 

gives the final result 

E+[l/ifn « (P„,(-A)<z_,)(0). 

□ 

When applied to the positive binomial distribution, this Theorem yields the following expansion (substi- 
tuting formula ([^0)1 for the factorial cumulants): 

Corollary 1. For K a positive binomial variate K ^ Bin(iV, p), the m-th Poisson approximation of the r-th 
inverse moment is given by 

m-l k / ly 

E+[l/K^] « g_.(0) + E ]^ E t^'^Hi-^y+'q-rm. (34) 

fe=l J=l ■'' 

Here, the coefficients aij obey the recurrence 

I 

"W=E7T^' (35) 

with boundary conditions q;o,o = 1 o.i^'d oiifl = for I > 0. 

Note that when fixing fi, this expansion is a series expansion in 



Proof. With = -N{j - (HZ]) becomes 

f N °° 1 \ 
t ^ expf--5]-(tpV)M 

oo oo 



j=o ■' 1=0 



Here we have defined the series expansion 



Et =E 

\k=2 J 1=0 



It is easily checked from this definition that the coefficients ai,j obey the recurrence stated in the corollary. 
The m-term approximation polynomial Pm is now obtained by imposing the constraint j + I < m — \ and 
setting t = 1, giving 

m— 1 , m—l—j 
]=0 1=0 

Substituting p — i-i/N and collecting terms in identical powers of N gives 

m-l , sj m-l-j 
3 = 1 ^' 1=0 

- i+EirE4^"-«(/-)^"'^- 

k=l 3 = 1 

To obtain the last line we have set k = j + 1 and reorganised the double summation. Combining this formula 
with ((33|) then gives the formula of the corollary. □ 

Table[T]gives the first values of aij. One sees that ais = + 2) and = 2[Hi+2 — + 4), where 
Hn is the n-th harmonic number. We are not aware of any closed form expression for j > 2. 
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Table 1: Values of the coefficients a/.j used in Corollary 1, for j + I < 7. 

For the special case r = 1, we can present an even more explicit formula. 

Corollary 2. For K a positive binomial variate K Bin(iV, p), the m-th Poisson approximation of its first 
inverse moment is given by 



m-l / k 

E+[i/x]«yo + E]^ E 

fc=l 



where 



1=1 



^^,^e-^Eli^,) + J2(^-ly.(e-^r)-l)^,-\ (37) 



and fi = Np. 

Proof. Equations and ([5^ yield explicit formulas for the sequence 

M-1,0 = e"^Er(^) 

M-i,a = |l-e 

Then this gives 

((-A)«g_i)(0) = e-^Er(/.)-^('M^^^^|-e-'^+^(-^)Vj! 





a-l 



a-l 



The third term simplifies, upon setting I = a — j and rearranging the double sum: 



a=l ^ / ^ j=Q a=l j=0 ^ ^ ■' 



1=1 j=0 



J + IJ j! 

I,,-' 



That yields 



((-A)>_i)(0) = e-^Er(A.) +E J - 1 

Substituting this expression for ((— A)"g_i)(0) in ([M)) of Corollary 1 gives the desired result. □ 



n\ ^\ {I - 1)! 



To illustrate the behaviour of the expansion, we depict the absolute and relative error of the approximation 
of E+[l/i4r] in Figures m and [5l respectively, for N = 10 and N = 100. It is clear from these figures that 
in contrast to previous expansions, the error is bounded uniformly over the complete interval [0,1]. The 
graphs have been produced using a Mathematica program, listed in Appendix [Cl We have also calculated 
the absolute error for the alternative series expansion p6p. It turned out that this expansion converged 
much more slowly than a Barbour expansion of the same degree (let alone one of the same order), even 
though the latter contains fewer terms than the former. For want of a better explanation, we attribute this 
phenomenon to the magic of Stein's method. 

The results depicted in Figure|4]can be compared to some explicit error bounds in |3]. Corollary 2.4 in [3] 
gives an upper bound to the absolute error when approximating expectations of a sum W oi N Bernoulli ('0- 
1') variates Xi. With pi = P[Xi = 1], the absolute error of the m-th order approximation to the expectation 
E[/i(W)] is bounded as 

^ i 

where /i = K[W] — J^iPi- ^ binomial variate K ^ Bm{N,p) is just a special case of this, obtained by taking 
all Pi equal. For the inverse moments, h is given by ft. = (0, 1, 2^'', 3^'', . . .), so that \\h\\ = 1. This gives the 
following bound: 

< 22"-ii^^7Vp™+i 

= 22™-i(l -e-^)p™. (38) 



Obviously, this bound is only useful for p < 1/4, and it matches the actual convergence only for p <ti 1- 
Nevertheless, this bound partially proves our claim that the Poisson approximation to the inverse moments 
converges. For bigger values of p, we currently have to rely on the numerical calculations reported in Figure 

SI 

A Numerical calculation of the inverse moments of a positive Pois- 
son variate 

In this Appendix we present a numerical method for calculating the inverse moments frin) = IE^[1/Q'^] of 
a positive Poisson variate with mean value fi. Series expansions are given in [Ml [TBI [27l [29l [32] but as these 
are asymptotic series they are not universally applicable. Plotting fr{p)/ reveals two different regimes in 
the range of ^i. For small p this function is seen to behave as for all r (see Figure[6|). This is also clear 
from the definition of /r, as the first terms of its defining series are 

p 2.2! ^ 3.3! 

This suggests that for small values of /i, the truncated series should give a good approximation: 



k 



fc=i 

For larger /i, the Figure suggests a 1/ ^"^ behaviour, with a moderately sharp cross-over region. For these 
larger values we can use the asymptotic series of |16) : 

M2-I I (r) I 



friP) - 2. T^' (40) 

t=0 ^ 



(k) 

where Sj are the Stirling numbers of the first kind. 

By appropriately choosing the cross-over point /i* at which to switch from ([39| to ([40)l and the number 
of terms Mi and M2 in the two series, one can tune the maximal relative error of the approximation while 
keeping the computational effort at bay. Tables [2] and [31 show the values for the cross-over point /j,*, and Mi 
and M2, as function of r, needed to obtain an approximation with relative error (|l-approximation/exact 
valuej) below 10~^ and 10~^°, respectively. In general, the choice of values for Afi and M2 involves a 
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Table 2: Values for the cross-over point /i*, and optimal number of terms Mi and M2, as function of r, to 
obtain an approximation of with relative error below 10^^. 



r 


1 


2 


3 


4 


5 


6 




25.734 


29.206 


33.99! 


B 37.903 


42.573 


47.068 


Ml 


63 


67 


74 


79 


85 


90 


M2 


20 


26 


33 


39 


46 


53 



Table 3: Same as Table[2]but for a relative error below 10 

trade-off between the two series. In the present case, however, the choice of M2 is determined because the 
second series is an asymptotic expansion. From a certain number of terms onwards, the benefit in including 
additional terms becomes marginal and ultimately the series diverges. We have chosen the value of M2 that 
minimises the value p^ below which the relative error becomes larger than the set minimum, so that the 
number of terms Mi of the first series, covering the remaining interval, can be made as small as possible. 



B Proof of the Poisson Expansion Formula 

In this Appendix we give a simple proof of the identity (jlSp underlying the Poisson expansion. Let / be a 
semi-infinite sequence 

/=(/(0),/(l),...,/(fc),...). 
As everywhere in this paper, we set f{k) ~ for < 0. Let / have mean value and let its fc-th factorial 
cumulant be with generating function 

( PC \ oo (j) 

fc=0 / j=0 

Recall that for any PDF /, = and k^^) = /i. 
Define the operator S 

M \ 



ki 

\k=2 



having matrix representation 



fc! 

\k=2 



Let g denote the sequence g = S'vr^ = S'.tt^, where tt^ is the sequence 

(71-^(0), 7r^(l), . . • ,7r^(fc), • . •) 
of the PDF of the Poisson distribution with mean value /i. We thus need to prove that g = /. 
The generating function of factorial cumulants can be written as 

log(f;/(fc)(l+x)M ^\ogfi{x), 

\fc=0 / 

where £,{x) is the semi- infinite vector 

for a; G C. This vector is the eigenvector of A corresponding to eigenvalue x. Applying A to the sequence 
^(x) indeed yields A^(a;) = (cc, x{l -\- x),x{l + x)^, ■ ■ ■) — x^{x). 

We will now calculate the factorial cumulant generating function of g. The inner product g'^£,{x) is given 

by 

Coo \ 

Since ^(x) is an eigenvector of A with eigenvalue x, this immediately gives 

Coo (fc-) \ 

The logarithm of the last factor is the factorial cumulant generating function of the Poisson distribution, 
which is known to be fix. Thus we get that the factorial cumulant generating function of g is 

log g^^ix) = ^^x'^+nx 

k=2 

- sV-'^ 

fe=0 

Since the right-hand side is identical to the factorial cumulant generating function of /, we have proven that 
/ = □ 



C A Mathematica program for the first inverse moment of a pos- 
itive binomial variate 



Here we reproduce the short Mathematica program that we have used to prepare Figure |3) For ease of 
implementation, the inverse moments of the Poisson distribution are calculated directly using Mathematica's 
summation capabilities, rather than via any recurrences or explicit formulas like the one of Proposition 1. 

(* First inverse moment of a-shifted Poisson: *) 

invmom [mu_ , a_] : = Sum [Exp [-mu] mu~k/k ! / (k+a) , {k , If [a==0 ,1,0] , Inf inity}] 

(* First inverse moment of 1-th forward difference of Poisson: *) 
invdif [mu_ , 1_] := Sum [Binomial [l,j] (-1) ~j invmom[mu, j] ,{j ,0,1}] 

(* Factorial cumulants of binomial: *) 
kappa[n_, p_, k_] = -n(k - l)!(-p)~k; 

(* m term Poisson expansion :*) 
expansion [n_ , mu_, x_ , m_] := 
Collect [Normal [ 

Series [Exp [Sum [kappa [n, mu/n, k] (-x t)~k/k!, {k, 2, m}] /t] , 

{t, 0, m - 1}]] /. t -> 1, x] ; 

(* The m-th order Poisson approximation (m=l,...,6) to 

the first inverse moment of Bin(n, p) :*) 
appr[n_, p_] = 

Table [(expansion [n, n p, x, m] /. x~k_ -> invdif [n p, k]) - 1 + 
invdif [n p, 0], {m, 1, 6}]; 

(* Exact expression: *) 

exact [n_ , p_] = Sum [Binomial [n, k]p~k(l - p) ~ (n - k)/k, {k, 1, n}] ; 
(* Absolute error: *) 

abserr[p_, n_, m_] := Abs[appr[n, p] [ [m] ] -exact [n, p] ] 
(* Relative error: *) 

relerr[p_, n_, m_] := Abs[l - appr[n, p] [ [m] ] /exact [n, p] ] 

(* Produces the graph of Fig. 4 for n=10, m = 1 to 6: *) 
<< Graphics ' Graphics ' 
grlist = Table [ 

LinearLogListPlot [ 

Table [{p, N [relerr [SetPrecision [p , 30], 10, m] , 30]} /. 

p -> k/500, {k, 1, 500}], PlotRange -> All, 
Plot Joined -> True], {m, 1, 6}]; 
Show[grlist, PlotRange -> {-12, 0}, DefaultFont -> {"Times-Italic", 16}, 
AxesLabel -> {"p", "rel . err . "}] ; 
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Figure 4: Absolute error as a function of p of the Poisson expansion of the first inverse moment of a positive 
binomial variate for A'' = 10 and A'' = 100, with 1 term (upper curve), and up to 6 terms (lowest curve). 
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Figure 5: Same as Figure 0] but showing the relative error. 
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Figure 6: Figure of fr{l^)/l^ = X^fcli ktw^ values of r from 1 (upper curve) to 6 (lower curve). The 
dashed curve represents the function e~^^. 



